Advection-diffusion

We first discuss how to solve

\begin{align*} u_t + \alpha u_{x} &= \kappa u_{xx}, \quad a < x < b, \quad t > 0,\\ u(a,t) &= g_0(t),\\ u_x(b,t) &= 0,\\ u(x,0) &= \eta(x). \end{align*}

One approach would be to combine everything and use one (implicit) method to time step the MOL discretization. Here one might want to use something like TR-BDF-2 (as opposed to trapezoid) so that if the eigenvalues from advection that lie on the imaginary axis leak into the right-half plane, after being perturbed, the method is still most likely stable. Although trapezoid is probably sufficient.

To illustrate splitting methods, we take a different approach and treat each term separately. In this lecture we won't do Strang splitting. We will just use the first-order (in time) approach (11.17).

Matrices for diffusion with Dirichlet at the left and Neumman at the right

Lax-Friedrichs and diffusion

Recall the setting of Lax-Friedrichs:

$$U^{n+1} = U^n - k A U^n + \frac{h^2}{2} B U^n$$

where $B$ is the second-order diffusion matrix. But we already have a matrix set up for diffusion in our problem so we could use that.

There is a problem though. As the advection parameter $\alpha$ increases, our required step size $k \leq h/|a|$ decreases relative to $h$. And as $k$ decreases relative to $h$, the effect of the diffusion term $h^2/2B$ will become more and more pronounced. We then have to decrease $h$ further to maintain the same accuracy. We have an eye toward going to two spatial dimensions and this will be way too costly. So, we need to use the Lax-Wendroff approach.

Lax-Wendroff diffusion and boundary conditions

Recall the derivation of the Lax-Wendroff method. If our MOL system with boundary conditions is given by

\begin{align*} U'(t) = A U(t) + g(t). \end{align*}

Here $A = - \frac{\alpha}{2h}S$ where $S$ is a matrix with $-1$'s on the first subdiagonal and $1$'s on the first superdiagonal. And $g(t) = \frac{\alpha}{2h}g_0(t) e_1$

To derive Lax-Wendroff with a Dirichlet condition at the left, we first ignore $g$

\begin{align*} U(t+k) &= U(t) + k U'(t) + \frac{k^2}{2} U''(t) + O(k^3)\\ &= U(t) + k A U(t) + \frac{k^2}{2} A U'(t) + O(k^3)\\ & = U(t) + k A U(t) + \frac{k^2\alpha^2}{8h^2} S^2 U(t) + O(k^3). \end{align*}

Then the key step in the Lax-Wendroff derivation is to replace $\frac{1}{4h^2}S^2$ with $B$, giving

\begin{align*} U^{n+1} = U^n + k A U^n + \frac{k^2 \alpha^2}{2} B U^n. \end{align*}

At a typical grid point this becomes

\begin{align*} U_j^{n+1} = U_j^n - \frac{\alpha k}{2 h} \left( U_{j+1}^n - U_{j-1}^n \right) + \frac{\alpha^2 k^2}{2 h^2} \left( U_{j-1}^n - 2 U_j^n + U_{j+1}^n \right). \end{align*}

Evaluating this at $j =1$ and using $U_0^{n} = g_0(t_n)$ tells us how to impose the boundary conditions.

We then can simulate a rudementary "wave tank" by decreasing the diffusion and turning off the boundary condition after some time.

Viscous Burgers' equation

We now discuss how to solve the nonlinear problem

\begin{align*} u_t + \alpha u u_{x} &= \kappa u_{xx}, \quad a < x < b, \quad t > 0,\\ u(a,t) &= g_0(t),\\ u_x(b,t) &= 0,\\ u(x,0) &= \eta(x). \end{align*}

Now, if we tried to combine everything into a single implicit method, we would have to solve a nonlinear system of equations at each time step. The splitting approach avoids this. Handling the diffusion is therefore the same. We have to think more carefully about how to handle the nonlinear advection term. We want a Lax-Wendroff-like approach but just replacing $\alpha$ with $\alpha u$ does not give a second-order method.

We follow the approach just below (10.18) in the text and just Taylor expand:

\begin{align*} u_t &= - \alpha u u_x,\\ u_{tt} & = - \alpha u_t u_x - \alpha u u_{xt}\\ & = 2\alpha^2 u u_x^2 + \alpha^2 u^2 u_{xx}\\ u(x,t+k) &= u(x,t) + k u_t (x,t) + \frac{k^2}{2} u_{tt}(x,t) + O(k^3),\\ u(x,t+k) &= u(x,t) - \alpha k u(x,t) u_x (x,t) + \frac{k^2}{2} \left( 2\alpha^2 u(x,t) (u_x(x,t))^2 + \alpha^2 (u(x,t))^2 u_{xx}(x,t)\right)+ O(k^3) \end{align*}

We now replace derivatives with centered differences:

\begin{align*} U^{n+1}_j = U^n_j - \frac{\alpha k}{2h} U_j^n \left( U_{j+1}^n - U_{j-1}^n \right) + \frac{\alpha^2 k^2}{4h^2} U_j^n \left( U_{j+1}^n - U_{j-1}^n \right)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_j)^2 \left( U_{j-1}^n - 2 U_j^n + U_{j+1}^n \right). \end{align*}

Setting $j =1$ gives the relation for the boundary:

\begin{align*} U^{n+1}_1 &= U^n_1 - \frac{\alpha k}{2h} U_1^n \left( U_{2}^n - g_0(t_n) \right) + \frac{\alpha^2 k^2}{4h^2} U_1^n \left(U_{2}^n - g_0(t_n) \right)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_1)^2 \left( g_0(t_n) - 2 U_1^n + U_{2}^n \right) \end{align*}

We write this as \begin{align*} U^{n+1}_1 &= F(U_j^n,U_{j+1}^n) + \frac{\alpha k}{2h} U_1^ng_0(t_n) - \frac{\alpha^2 k^2}{2h^2} U_1^n U_2^n g_0(t_n) + \frac{\alpha^2 k^2}{4 h^2}U_1^ng_0(t_n)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_1)^2 g_0(t_n). \end{align*} or in matrix-vector form: \begin{align*} U^{n+1} &= U^n - k\,\mathrm{diag}(U^n) (A U^n) + k^2 \mathrm{diag}(U^n) (A U^n)^2 + \frac{\alpha^2 k^2}{2}\mathrm{diag}(U^n)^2 (B U^n) + h^n e_1,\\ h^n &= \frac{\alpha k}{2h} U_1^ng_0(t_n) - \frac{\alpha^2 k^2}{2h^2} U_1^n U_2^n g_0(t_n) + \frac{\alpha^2 k^2}{4 h^2}U_1^ng_0(t_n)^2 + \frac{\alpha^2 k^2}{2 h^2} (U^n_1)^2 g_0(t_n). \end{align*}

Neumann BCs

Waves in $(2 + 1)$-dimensions